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ABSTRACT 

This  thesis  examines  the  differences  of  deterministic  and  stochastic 

LANCHESTER-type  combat  models.  Using  an  example  of  square-law  attrition, 
solution  methods  and  solutions  are  described.  A  new  analytic  solution 
for  equal  attrition  rate  coefficients  is  given.  The  numerical  comparison 
includes  hypotheses  about  the  expected  force  levels  and  the  variability 

in  the  expected  force  levels  as  a  function  of  time,  initial  force  levels, 
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I.   INTRODUCTION 

Combat  models  are  widely  used  as  decision  aids  in  the  defense- 
planning  process,  at  least  within  the  NATO  alliance.  Current  operation- 
al combat  models  are  wery   complex  because  combat  is  a  very  complex 
process.  Unfortunately  it  is  difficult  (if  not  impossible)  for  the 
beginner  to  understand  the  modelling  approaches,  concepts  and  motivation, 
that  may  have  been  used  to  build  such  operational  models.  However,  one 
frequently  considers  a  simple  model  as  a  paradigm  for  the  development  and 
understanding  of  such  complex  models.  This  basic  approach  will  be  used 
in  this  thesis  to  explore  certain  issues  in  the  on-going  debate  about 
the  relative  merits  of  stochastic  and  deterministic  combat  models. 

A  simple  model  is  examined  to  explore  differences  between  a  determin- 
istic and  a  stochastic  approach  to  a  certain  type  of  analytical  combat 
model.  As  already  mentioned,  combat  is  a  very  complex  process,  but  it 
is  also  a  complex  random  process,  which  can  be  supported  by  many  examples 
from  military  history.  Analytical  models  are  abstractions  and  yery   often 
simplifications  of  reality.  It  seems  to  be  a  legitimate  question  to  ask, 
what  effects  the  further  abstraction  of  neglecting  the  randomness  in 
combat  may  have.  At  this  moment,  it  should  be  pointed  out  that  within 
existing  operational  analytical  models,  both  stochastic  and  deterministic 
models  are  used. 

Previous  work  done  by  SPRINGALL  [9]  and  CLARK  [4]  evolved  around 
theoretical  aspects.  Their  main  concern  was  to  give  exact  analytical 
solutions  and  their  proofs.  CRAIG  [5]  started  to  explore  the  differences 


between  stochastic  and  deterministic  models  more  from  the  numerical 
point  of  view,  which  will  be  continued  in  this  thesis. 

In  the  next  chapter,  a  deterministic  and  stochastic  version  of  a 
differential  combat  model  will  be  described.  The  deterministic  versions 
are  well-known  as  LANCHESTER's  equations  of  modern  warfare,  which  were 
developed  in  1914.  Combat  models,  which  model  attrition  from  enemy 
action  through  a  system  of  differential  equations,  are  usually  referred 
to  as  LANCHESTER-type  models  of  warfare. 


II.  THE  PARADIGM  MODELS 

A.  THE  DETERMINISTIC  MODEL 

First,  LANCHESTER's  equations  of  modern  warfare  will  be  briefly 
reviewed  and  some  simple  extensions  given. 

In  1914  LANCHESTER  [7]  hypothesized  that  under  "modern  conditions" 
in  a  combat  between  two  homogeneous  forces  the  firepower  of  the  surviving 
weapons  of  one  side  can  be  concentrated  on  the  surviving  targets  of  the 
other  side,  so  that  each  side's  casualty  rate  is  proportional  to  the  num- 
ber of  enemy  firers.  This  can  be  described  by  the  following  equations: 

^=-ay  (2.1) 

^=-bx  (2.2) 

with  initial  conditions 

X(0)  =  Xq  (2.3) 

Y(0)  =  yQ  (2.4) 

where  a  is  the  attrition  rate  with  which  the  Y-force  attrits  the 
X-force,  similarly  for  b.  Xq  and  Y^  are  the  initial  force  levels  and 
X(t)  and  Y(t)  are  the  force  levels  at  time  t.  The  force  levels,  as  a 
function  of  time  t,  can  be  written  as 

X(t)=XQCOsh(Vab  t)-Va7b  yQSinh(Vab'  t)  (2.5) 

Y(t)=yQCOsh(^ab'  t)-Vb7a  XQSinh(Vab  t)  (2.6) 

The  state  equation  relating  initial  force  levels  with  force  levels  at 
some  time  t  can  be  derived  by  dividing  (2.1)  by  (2.2),  separating 
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variables  and  integrating  to  yield 

b(xQ2-X(t)^)=a(yQ2.Y(t)2).  (2.7) 

This  form  of  the  state  equation  explains  why  this  model  is  also  referred 
to  as  the  "square-law"  attrition  process.  WEISS  [111  has  given  a  set  of 
assumptions  under  which  LANCHESTER's  equations  for  modern  warfare  may 
apply: 

Al . )  Two  homogeneous  forces  are  engaged  in  combat.  Every  unit 
on  a  particular  side  has  the  same  capabilities.  The  attrition 
rate  may  be  different  for  the  two  forces. 

A2.)  Each  unit  on  one  side  is  within  weapon  range  of  all  units 
on  the  other  side. 

A3.)  The  effects  of  successive  rounds  on  the  target  are 
independent. 

A4.)  Each  unit  has  perfect  knowledge  of  target  locations  and 
fires  only  at  live  target  (one  at  a  time)  killing  them  at  a 
constant  rate,  which  does  not  depend  on  the  number  of  targets 
alive. 

A5.)  Fire  is  uniformly  distributed  over  surviving  targets. 

The  above  model  implies  a  fight  until  one  force  is  annihilated. 
Therefore  the  model  will  be  slightly  changed  by  introducing  the  concept 
of  unit  breakpoints,  X.   and  Y.  ,  which  are  simply  force  levels  at  which 
the  side,  who  reaches  it  first  "breaks  off"  the  engagement,  leaving  the 
other  side  as  a  winner.  Also,  to  be  more  precise,  it  should  be  noted 
that  negative  force  levels  for  breakpoints  equal  zero  or  force  levels 


less  than  nonzero  breakpoints  are  impossible.  So  the  deterministic 
LANCHESTER-type  combat  model  with  "square-law"  attrition  takes  the 
following  form: 


dx  _  . 
Ht"  ^ 


dt 


-ay  X|^p<X(t) 

0    otherwise 
-bx  Xj^p<X(t) 

ybp<Y(^) 

0    otherwise 


(2.8) 


(2.9) 


with  initial  conditions 

X(0)  =  Xq  (2.10) 

Y(0)  =  yQ  (2.11) 

The  model  in  this  form,  equations  (2.8)  through  (2.11),  was  used  for 

comparisons  throughout  the  thesis. 

B.  THE  STOCHASTIC  MODEL 

There  are  several  ways  to  include  random  variations  in  LANCHESTER- 
type  models.  These  are: 

*  The  attrition  rate  coefficients  may  be  random  variables. 

*  The  enemy's  initial  force  level  may  be  a  random  variable, 
weakening  the  assumption  of  perfect  knowledge. 

*  The  breakpoints  may  be  random  variables. 

*  The  casualty  rate  is  fixed,  but  the  occurrence  of  casualties 
over  time  may  be  random. 

The  only  random  variation  considered  here  will  be  the  random  occur- 
rence of  casualties  over  time.  Another  specification  was  to  choose  a 
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model  similar  to  the  "square-law"attrition  in  order  to  allow  comparisons. 
In  other  words,  the  question  to  be  asked  is  "how  do  random  fluctuations 
in  the  occurrence  of  casualties  modify  the  deterministic  results  of  the 
square-law  attrition  process?" 

The  approach  used  here  was  a  continuous  parameter  MARKOV  chain  model, 
where  the  time  t  varies  continuously  and  the  number  of  combatants  on 
each  side  is  a  non-negative  integer.  Let  M(t)  be  the  size  of  the  X-force 
at  time  t  with  a  particular  state  value  m.  Let  N(t)  be  the  size  of  the 
Y-force  at  time  t  with  a  particular  state  value  n.  Let  m^  and  nQ  be  the 
initial  force  levels  and  m.  ,  n.   be  the  breakpoint  force  levels  of  X 
and  Y  respectively.  Fig.  1  shows  the  state  space  of  this  MARKOV  chain 
model.  Note  that  at  a  given  time  t,  any  state  is  described  by  the  two 
force  levels  of  the  X  and  Y  force.  As  each  side  loses  units  due  to 
attrition  and  no  replacements  are  allowed,  it  is  easy  to  understand  why 
BILLARD  [11  referred  to  this  type  of  process  as  a  "bivariate  death  process." 


n 


(."^o,%) 


Figure  1  -  STATE  SPACE  OF  THE  MARKOV  CHAIN  MODEL 
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For  the  description  of  the  stochastic  square-law  attrition  process 
corresponding  to  the  two  deterministic  differential  equations  (2.8)  and 
(2.9),  a  system  of  many  differential  equations,  depending  on  the  battle 
termination  model,  is  required.  This  system  will  be  given  for  a  fixed- 
force-level  -breakpoint  battle  with  square-law  attrition. 

The  following  assumptions  yield  the  stochastic  square-law  attrition 
process. 

Al . )  The  attrition  process  depends  only  on  the  current  system 

state  and  time  and  not  on  the  past  history  (this  assumption  is 

usually  referred  to  as  markovian  property). 

flo  ^  Tu«  ov,oK=k-ii-;+w  of  one  x  casualty  during  the]  _  . 
A2.)  The  probability  P^^.^^  interval  t  to  t+h   i    "  ^^ 

A3  )     ThP  Drobahilitv  P  f °"^  ^  casualty  during  the  \  =  uu 
M.)     ihe  probability  i^^^ime  interval  t  to  t+h   /   °" 

A4.)  The  probability  of  more  than  one  casualty  occurring  in 
the  time  interval  t  to  t+h  is  of  the  order  of  magnitude  o(h), 
where  limu^o(h)/h  =  0. 

A5.)  No  more  casualties  can  occur  once  m  =  m.  or  n  =  n.  . 

Making  the  time  interval  h  infinitesimally  small,  the  following 
set  of  forward  CHAPMAN- KOLMOGOROV  equations  can  be  developed.  Let 
P(t,m,n)  be  the  probability  that  the  system  is  in  state  (m,n)  at  a 
time  t.  For  convenience  each  equation  is  related  to  a  region  in  the 
state  space  shown  in  Fig.  2. 
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For   m=mp,    and    n  =  n«,    Region    I 


"0 


0 


dP 

^(t,mQ,nQ)    =    -(anQ+brriQ)    P(t,mQ,nQ)  (2.12) 

for  "ikr,<"^<'^n   ^"^   "~"n'    ^^Qio"    II 

^(t,m,nQ)    =    anQ    P( t ,m+l .ng)- (anQ+bm)    P(t,m,nQ)       (2.13) 
for   nun<n<nQ   cind   m=mQ ,    Region    III 

dP 

j^(t,mQ,n)    =    bniQ    P(  t  ,mQ  ,n  +  l  )- (an+bmQ)    P(t,mQn)         (2.14) 

for  m.  „<m<m-.    and    n.  „<n<n«,    Region    IV 
bg 0 _bjg 0  ^ 

dP 

^(t,m,n)    =    anP(t,m+l ,n)+bmP(t ,m,n+l )-(an+bm)P(t,m,n) 

(2.15) 

for  m=m.  „   and    n.  „<n<n«.    Region    V 
D_2 b_g 0  ^ 


d  P 

jf(t,m^p,n)    =    anP(t,mj^p+l  ,n) 

for   n  =  n.  „    and   m. „<m<mn.    Region    VI 
b_2 b_2 0  ^ 

dP 

^(t,m,njjp)    =    bmP(t,m,n|^p+l  ) 

for  m=m.  „   and    n  =  n. „,    Region    VII 
D_2 bp  ^ 

P(t,m^p,njjp)  =  0  for  all  t 


(2.16) 


(2.17) 


(2.18) 


because  of  the  definition  of  a  breakpoint  force  level.   The 
initial  condition  is 


P(0,mQ,nQ)  =  1  .0 


(2.19) 
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Figure  2  -  REGIONS  IN  THE  STATE  SPACE 
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As  P(t,m,n)  is  a  joint  probability  distribution,  the  follow- 
ing must  also  be  true 


0^  P(t,m,n)^l .0 


(2.20) 


and 


itIq     nQ 
bp     bp 


P(t,m,n)  =  1 .0 


(2.21) 
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III.   SOLVING  THE  DIFFERENTIAL  EQUATIONS 

A.  SOLUTIONS  FOR  THE  DETERMINISTIC  MODEL 

Force  levels  as  a  function  of  time  were  already  given  in  equations 
(2.5)  and  (2.6).  It  is  relatively  easy  to  obtain  analytical  solutions 
for  such  simple  deterministic  models  as  described  before.  On  the  other 
hand,  it  is  necessary  to  point  out  that  for  models  with  any  degree  of 
operational  realism,  analytical  methods  for  solving  the  differential 
equations  are  usually  not  available.  Therefore,  some  numerical  method 
with  a  digital  computer  is  usually  used.  TAYLOR  [10]  has  summarized  in 
his  Appendix  C  the  most  widely  used  numerical  methods,  a  discussion  of 
which  seems  unnecessary  at  this  point. 

B.  SOLUTIONS  TO  THE  STOCHASTIC  MODEL 

Even  for  this  relatively  simple  stochastic  model  with  fixed  force 
level  breakpoints,  which  are  usually  nonzero,  a  complete  set  of  general 
analytical  solutions  for  the  differential  equations  (2.12)  through  (2.19) 
has  not  been  found.  On  the  other  side,  given  some  minor  restrictions 
like  a  fight  to  the  finish  or  equal  attrition  rate  coefficients,  solu- 
tions, or  at  least  solution  methods  have  been  proposed  which  will  be 
briefly  discussed  in  the  next  section. 

First,  the  method  for  getting  the  state  probabilities  used  here  will 
be  described.  Numerical  solutions  were  obtained  using  the  fourth-order 
RUNGE-KUTTA  method,  which  is  probably  one  of  the  best  known  finite  dif- 
ference approximations  to  ordinary  differential  equations  (next  to  the 
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EULER-CAUCHY-method).  To  increase  the  accuracy  of  the  overall  solutions 
analytical  results  for  certain  regions  of  the  state  space  were  substituted, 
These  analytical  solutions  will  be  stated  now.  For  region  I,  i.e.  no 
casualties  on  either  side,  the  solution  to  (2.12)  is 

P(t,mQ,nQ)  =  exp  -(anQ+bniQ)t  (3.1) 

which  can  be  derived  by  the  standard  method  used  for  this  kind  of  dif- 
ferential equation.     For  the  boundary  cases,  region  II  and  III,  i.e. 
one  of  the  two  sides  has  not  yet  had  a  casualty,  TAYLOR   [101   has  given 
the  analytical  expression  as 
for  nikn^^^^'^O    ^"^   "~"0'    f^^gion    II 

P(t,m,nQ)    =   i    |anQ/b(e'^^-l)l^    exp(  -  (biriQ  +  anQ)    t)    ,      (3.2) 
where   J    =   niQ-m 
for   nkn<n<nQ    cind   m=mQ,    Region    III 

P(t,mQ,n)  =  1  |bmQ/a(e^^-l)l^  exp(- (brnQ+ang)  t)  ,   (3.3) 

where  K  =  nQ-n. 

These  two  equations  (3.2)  and  (3.3)  were  obtained  by 
recursively  solving  equations  (2.12),  (2.13)  and  (2.14)  "from 
the  top  down . " 

C.   OTHER  ANALYTICAL  SOLUTIONS 

The  solutions  or  solution  methods  for  getting  the  state 
probabilities  will  only  be  stated  for  the  square-law  attri- 
tion process.   Only  two  were  used  for  the  numerical  work  for 
this  thesis. 
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Apparently,  one  of  the  "oldest"  analytic  solutions  was 
given  by  BROWN  [3]  in  1955  for  the  general  stochastic 
LANCHESTER-type  combat  model  with  time  independent  attrition 
rates.   His  approach  and  solution  will  be  briefly  outlined 
for  square-law  attrition.   Consider  a  path  from  state 
(  mQ  ,  n^  )  to  some  state  (m,n).   This  path  can  be  described 
as  a  sequence  of  J  =  (  m^  -m)  zeros  and  K  =  (  nQ  -n)  ones, 
where  a  zero  corresponds  to  an  X  casualty  and  a  one  to  a  Y 
casualty.   Using  the  binary  representation  of  a  positive 
integer,  one  can  define  to  each  realization  of  a  battle  path 
an  integer  k  given  by 

^   "  ^k,l^l<,2-  •  •^k,J  +  K 

a.  L. 

where  d.    =  1  if  the  r   casualty  along  a  battle  path 
K ,  r 

corresponding  to  k  is  a  Y  casualty  and  cl.    =  0  otherwise. 
Also  let  Ij  |.  be  the  set  of  all  positive  integers  whose 
binary  representation  contains  exactly  K  ones  and  J  zeros. 
Then 


,-.•£ 


""k.r  =  '"o"'-^  fr^  \,i 


r 


"k,r  =  "O   -  -^  ^k,j 


Then  BROWN  [31  has  shown  that 
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-  J+K-1 


kel      '^         -  -^ 


k.r 


J,K 


r=0 


exp(-iut)-exp(-tl  (m|^^^,n,^^^)) 


du 


,  (3.4) 


where  i  =  y-l  , 

1 (m  ,n )  =  an+bm  , 

lk,r+l 


K 


k,r   1-iu/  Um|^^^,n,^^^) 


and 

9k,r+l  =  ^k,r+l^"k,r'-^^-^k,r+l^'^"^k,r  * 
There  was  no  indication  that  nonzero  breakpoints  were 
excluded.   A  discussion  of  this  solution  follows  in  the  next 
section . 

About  14  years  later,  in  1969,  CLARK  [4]  proposed  another 
approach  which  TAYLOR  [10]  called  a  "hybrid  analyti cal -numer- 
ical" method.   The  restriction  is  that  the  breakpoints  have 
to  be  zero,  i.e.  it  is  a  fight  to  the  finish.   Although 
proposed  for  a  general  time  independent  attrition  function, 
this  approach  will  be  outlined  for  "square-law"  attrition. 
Then  according  to  CLARK  [4]  the  state  probabilities  are  given 
by 
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for   0<m<m^   and   0<ninQ 


m„        n 


"'0        "iL. 
P(t,m,n)    =       Y^     Y2    ^j!k    •exp[-(al<  +  bj)t)  ,    (3.5) 


j=m     k=n 


for   0<m<mj^   and 


P(t,m,0)    = 


n  =  0 


niQ        nQ 


cS'o   ^21      21    C'l^l    .exp(-(ak  +  bj)t),    (3.6) 
^*^        j=m      k  =  l      ^'"^  ^  ^ 


for  m=0   and   0<n<n 


0 


P(t,0,n)    =    cg;^   + 


riiQ      nQ 


j=m   k=l 


C^l^    •exp(-(ak+bj)t)  ,(3.7) 


and  at  last  thereis  to  remember  that  P(t,0,0)  =  0  for  all 
times  t.   The  constants  c'J'J  are  determined  by  a  system  of 
partial  difference  equations. 

For  0<m<j<mQ  and  0<n<k<nQ 


.m ,  n 
'j,k 


-m+l »n   ^     pm  ,n  +  l 
an   j,k       bm  j,k 

a(n-k)+b(m-j) 


,  (3.8) 
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for  0<m<j<mp|  and  0<n  =  k<nQ 


m,n       m+1  ,n 
C  .   =  an  C  . 
J ,n       J  ,n 


b(ni-j) 


.    (3.9) 


for   0<ni  =  j<niQ   and   0<n<k<nQ 


m,n  m,n+l 

C„    .    =    bm   C 
m ,  k  m  ,  k 


a(n-k) 


,      (3.10) 


for  0<m=j<mQ  and  0<n=k<n, 


but    (m,n)    f    (mQ.n^) 


m,n  '"O        "'O  m,n         ,^IL-         ni,n 

C       =   -  2Z     211    c.    .     -  Z_      C.    . 
"•'"  j=m      k  =  n  +  l    ^'^  j=m+l      ^ '^ 


,      (3.11) 


with    ,"|0'"0    __    ^    0^ 
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Also  for  0<m<j^mQ  and  0=n<k<nQ 


Cjik  =   -  ^"^  J>k  .        (3.12) 


ak+bj 
similarly  for  0=m<j<mQ  and  0<n<k<n^ 


rO,n         pl,n 
J  »K         J  ,k 


ak+bj  .        (3.13) 


Then    for   l<m<mQ 

j=m      k=l 
and   finally    for   l<n5nQ 

m,0 
^j,k 

rriQ        nQ 

^S:S  =  ->    > 

j=l      k=n 

pO,n 

(3.14) 


(3.15) 


Though  having  the  publishing  date  of  1979,  the  next 
approach  was  published  in  June  1980  by  BILLARD  [1].   She 
considered  the  LANCHESTER-type  square-law  attrition  combat 
model  as  a  pure  death-process  and  applied  SEVERO's  [81 
recursive  theorem  for  solving  differential  equations.   As 
before,  only  a  fight  to  the  finish  has  been  considered. 
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The  first  step  is  to  identify  each  point  (m,n)  in  the 
state  space  by  a  counting  coordinate  k ,  where 


k  =  (  mQ+l)(  nQ+1) 
-m(  nQ+l)-n 


(3.16) 


Then 


P[t,m,n]=  P[t,k]  (3.17) 

and  the  differential  equations  (2.12)  through  (2.18)  take  on 

a  slightly  different  form.   As  an  example,  (2.15)  will  be 
given  by 

|^(t,k)  =  anP[t,k-nQ-ll+bmP[t,k-l] 

-(an+bm)P[t,k]   .        (3.18) 
The  whole  set  of  differential  equations  was  then 
expressed  in  matrix  terms  as 
d 


^P(t)  =  BP(t) 

with  a  solution  given  as 
P(t)  =  Ce  (t) 


(3.19) 


(3.20) 


where  ^(t)  is  the  (m^  +  l  )  (n^^+l  )xl  -vector  with  elements 
exp(b. t)  with  b.  being  the  k   diagonal  element  of  the 
matrix  B^.   The  matrix  B^  can  be  partitioned  into  submatrices, 
whose  m-coordinate  is  common,  due  to  the  ordering  defined  by 
the  counting  coordinate  k  (Equation  3.16). 

Then 

i  =  (^uv^'  "'^  "  "•  »2...mQ  +  l 


24 


where  the  submatrices  b    have  the  elements 
-uv  "  (b^^(p,q)),  p,q  =  l,2...nQ+l. 

So  the  matrix  B^  has  the  elements 
bM.,(P»P)  =  -a(n«-p+l  )-b(mn-u  +  l  ) 


for  u  =  l,2...mQ  and  p  =  l,2...n 


0 


b^^jj(p,p-l)  =  b(mQ-u+l)  for  u  =  l,2...mQ  +  l 


and 


p  =  2,3.  . .nQ  +  1 


^u,u-l^P'P^  "  a(nQ-p+l)  for  u  =  2,3...mQ+l 


P  =  1  ,2. .  .nQ  +  1 . 

All  other  elements  are  zero. 

Thus,  the  matrix  B^  has  at  most  3  nonzero  entries  per  row 

or  column.   The  matrix  C_  can  be  partitioned  in  the  same  way. 

Then  using  SEVERO's  [8]  theorem  and  the  special  form  of  the 

matrix  B^,  only  a  part  of  the  C_-matrix  needs  to  be  determined. 

This  part  will  be  omitted  here,  but  the  final  result  will 

be  given  by 

k 
P(t,k)  =  ^  c(k,j).exp(bjt)  (3.21) 

j  =  l 

J.  U 

where  c(k,j)  is  the  (k,j)    element  of  the  solution  matrix  C^. 
The  previous  two  approaches  have  required  that  the  force 
level  breakpoints  be  zero.   Now,  a  result  will  be  given  whose 
restriction  is  that  the  attrition  rate  coefficients  be  equal, 
but  nonzero  breakpoints  are  allowed.   For  further  reference 
it  will  be  called  the  Equal -Attri ti on-Rate-Coeffi ci ent- 
Solution  (EARCS). 
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Let  a  =  b  =  f. 


^or   m.  ^<m<m^  and  n. „<n^nn 
bp 0 b_g 0_ 


Pft  m  n^  =    C(m,n) 


(l.e-^t)"^0''"0-'"-" 


exp(-f (m+n) t) 


(3.22) 


for  m  =  m.  ^  and  n.  „<n<nn 
bp      bp    0 


P(t,m.  ,n)  =   ^"    C(m^p^1,n)  ^ 


bp 


±_   (-1)^ 
k  =  0 


J! 
/j\   Jl-exp(-ft(ni^p^Un^-k))] 


(2.23) 


where  J  =  m^+nn-m.  „-l -n 
0   U   bp 

for   n    =    n.  „    and   m.     <m<mp, 
bp  bp  0 


D/+   m   r,      ^    -      "^m      C(m  ,n.  „+l  ) 
P(t,m,n.     )    =      ^        bp      ' 

°P  K! 


(-1) 


j  =  0 


'l-exp(-ft(m+Un^p+j)) 
f(m+l+nj^p+j) 


}• 


(3.24) 


where    K   =   m^  +  n^-m-n.  „-l  . 
0      0  bp 


The    coefficients    C(m,n)    satisfy   for  mbrj<m<nio    cind    "hn'^^^^O 
the   following    partial    difference    equation 

C(m,n)    =    nC(m+l,n)    +   mC(m,n+l)  (3.25) 
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with    the   boundary   conditions 

C(m,nQ)    =    (nQ)  for   mj^p<m<mQ 


and 


C(mQ,n)  =  (ftIq)     for  nj^p<nlnQ 


This  result  has  been  developed  using  a  method  verbally  pro- 
posed by  TAYLOR.   The  method  will  now  be  outlined.   ISBELL 
and  MARLOW  [6]  described  a  stochastic  LANCHESTER-type  attri- 
tion process  with  a  different  attrition  function.   Instead  of 
the  attrition  of  one  force  being  proportional  to  the  number  of 
of  enemies  of  the  opposing  force  as  in  the  square-law  attri- 
tion (e.g.  for  the  attrition  of  the  M-force  let 

A(m ,n )  =  an 
be  the  attrition  rate  and 

B(m,n )  =  bm 
the  attrition  rate  for  the  N-force  with  square-law  attrition), 
their  attrition  rates  looked  like 

A(m,n)  =  an+cm 
and 

B(m,n)  =  bm+dn 
with  the  restriction  that 

a+c  =  b+d. 
But  with  c=d=0  and  a=b=f  we  are  back  to  square-law  attrition. 
This  leads  to  equations  (3.22)  and  (3.25).   Equations  (3.23) 
and  (3.24)  were  derived  in  the  following  way  (e.g.  for  3.23). 
Solving  equation  (3.22)  for  m  =  m^^  +1  yields 
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C(m.  „  +  l  ,n)         /  \ 


'0  "0  '"bp 


-  f  t .  ni  n + n  «  -  m .  „  - 1  -  n 
.  (i-e  ^^)  0   0   bp 


(3.26) 


substituting  for  the  second  factor  its  BINOMIAL  expansion 


k=0       ^  ' 


^-ftk(,)J-k 


with  J  =  "iQ'*'nQ-m.  -1  and  multiplying  through  by  the  third 
factor.   Then  using  equation  (2.16),  the  differential  equation 
for  m  =  m.   and  "kn'^^-^o  »  ^""^  substituting  equation  (3.26) 
into  the  extended  form,  it  can  now  easily  be  integrated  to 
yield  equation  (3.23). 

D.   DISCUSSION 

In  the  discussion  of  the  analytical  solutions  outlined  in 
the  last  section,  there  is  one  important  point.   BROWN  [3] 
himself  points  out  that  unless  m  is  close  to  mg  and  n  is 
close  to  nQ  his  result  (equation  (3.4))  is  of  "little  practi- 
cal interest."   Most  of  the  analytical  sol  utions ,  especi al ly  for 
more  general  LANCHESTER-type  models,  have  little  more  than 
"symbolic"  character.   BROWN'S  solution  is  a  good  example 
of  that. 

In  comparing  the  solutions  given  by  BILLARD  [2]  and 
CLARK  [4],  this  author  has  the  feeling  that  both  solutions  are 
equivalent  and  only  the  representation  is  different.   This 
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intuitive  guess  needs  verification.  It  may  be  concidence 
that  CLARK  [4]  and  SEVERO  [81  published  their  work  in  the 
same  year. 

The  last  presented  solution  (equation  (3.22)  through 
(3.24))  seems  to  be  relatively  handy  for  use  on  a  computer. 
It  has  a  big  advantage  over  numerical  solution  methods  other 
than  its  accuracy,  because  it  is  an  exact  result.   Like 
CLARK'S  method,  the  coefficients  have  to  be  calculated  only 
once  for  a  given  set  of  input  data.   Then,  to  get  the  state 
probabilities  for  a  certain  point  in  time  you  have  to  make 
only  one  set  of  calculations,  as  opposed  to  the  numerical 
methods  where  one  has  to  go  from  time  t  =  0  to  time  t  =  t 
in  small  time  steps  and  then  have  only  an  approximate  result 

E.   IMPORTANCE  OF  THE  STATE  PROBABILITIES 

The  state  probabilities  as  a  function  of  time  are  the 
key  to  calculating  several  quantities  of  interest.   These 
are  expected  force  levels  as  a  function  of  time,  variances 
and  standard  deviations  in  the  force  levels  and  also  the 
probability  of  winning.   These  quantities  are  necessary  to 
legitimately  compare  the  stochastic  with  the  deterministic 
results . 

To  get  at  least  a  feeling  of  how  the  state  probabilities 
evolve  over  time,  the  joint  probability  distribution  will  be 
presented  in  a  3-D-picture.   It  is  indeed  surprising  that 
more  use  has  not  been  made  of  computer  graphics  to 
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investigate  the  dynamics  of  a  stochastic  LANCHESTER-type 
combat  model.   Table  1  gives  the  data  used  for  the  next  five 
figures.   These  figures  may  be  thought  of  as  "snapshots"  of 
the  joint  probability  for  the  survivors  in  this  battle  taken 
at  a  sequence  of  increasing  times. 
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TABLE  1 


Data  for  the  Numerical  Example  1 


Force  Levels 


%  =  ^0 


"0  =  ^0 


bp 
bp 


Attrition  Rates 


a  =  0.008 


b  =  0.004 


M  casualties  per 
minute  and  N  firer 

N  casualties  per 
minute  and  M  firer 


At  Times  t^  =  0.025  t^ 

t2  =  0.250  t^ 
t3  =  0.500  t^ 
t^  =  0.750  t^ 
tg  =  1 .000  t^ 

where  t^  =  155.81  minutes  is  the  time  a 
deterministic  battle  with  the  same  force 
levels  and  unit  breakpoints  ends,  i.e. 
x(t^)  =  0.0  and  y(t^)  =  28.28 
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i  P(t,in,n) 


Pore*  L€v«l,  m 


Force  Level,  n 


Figure  5  -  FLOT  OF  JOINT  PROBABILITIES  P(t,M,n; 
with  data  according  to  Table  1 
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a  p(t,m,n) 


Fore*  L«v«l,    m 


Forca  L«vel,    n 


Jfigur*  if  -   fLOT  of  joint   PROBABILITIES   P(t,»,n; 
with  data  according  to   rabl*   1 
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Hp(t,s,n) 


Force  L«v«l,   m 


Force  Level,   n 


Figure   5  -  -PLOT  OF  JOINT   PROBABILITIES   P(t,«,n; 
with   data  according  to  Tabl«   1 
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Force  L«vttl,  n 


Force  Level,  n 


rigure  6  -  fLOT  OF  JOINT  PROBABILITIES  P(t,M,n) 
with  data  according  to  Tabl«  1 
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,jP(t,in,n) 


rorc*  Level,   m 


Force   Level,    n 


Jrlgur*  7   -   PLOT   OF  JOINT   PROBABILITIES   P(t,M,n) 
vlth   data   according  to  Tabl«   1 
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At  the  time  t  =  0  all  probability  is  located  at 
(mQ,nQ)  in  the  state  space  (Region  I).   As  time  passes,  the 
probability  mass  is  distributed  over  more  states,  with  the 
mode  moving  away  from  the  diagonal  towards  the  winning  side. 
All  points  in  the  state  space  with  breakpoints,  i.e.  (mkn>n) 
and  (m,n.  )  for  all  m  and  n,  are  absorbing  states,  probability 
mass  is  absorbed  in  that  states.   The  sum  of  probability  mass 
in  Region  V  (see  Fig.  2)  represents  the  probability,  that  the 
N-force  wins  at  that  given  time,  in  Region  VI  that  the  M-force 
wins. 

The  next  five  figures  show  a  similar  sequence  of  plots  for 
the  joint  probability  P(t,m,n)  for  the  force  levels  M(t)  and 
N(t).   The  data  is  explained  in  Table  2.   Note  the  small  dif- 
ferences because  of  the  nonzero  force  level  breakpoints. 
Probability  mass  having  reached  the  breakpoint  "piles"  up 
there.   The  state  space  is  reduced  by  the  fixed  force  level 
breakpoints,  but  the  probability  distribution  evolves  basical- 
ly in  the  same  qualitative  manner  as  in  the  previous  example. 
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TABLE  2 


Data  for  the  Numerical  Example  2 


Force  Levels 


mQ  =  40 

"0  -  ^0 


bp 
bp 


Attrition  Rates 


a  =  0.008 
b  =  0.004 


M  casualties  per 

minute  and  N  firer 

N  casualties  per 

minute  and  M  firer 


At  Times  t^  =  0.025  t^ 

t2  =  0.250  t^ 
t3  =  0.500  t^ 
t^  =  0.750  t^ 
tg  =  1 .000  t^ 

where  t^  =  120.68  minutes  is  the  time  a 
deterministic  battle  with  the  same  force 
levels  and  unit  breakpoints  ends,  i.e. 
x(t^)  =  8.0  and  y(t^)  =  28.83 
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p(t.in.n) 


force  Lev«l#    «» 


fore*  L«vel.    n 


,        ,U>.Or.OmPHOBABlUTXKSP(....») 
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i>  P(t,ni,n) 


Toxem  Level,    m 


Force   Level,    n 


Flgur«   9   -   PLOT   OF   JOINT   PROBABILITIES   P(t,B,n) 
with   d«ta   according   to    rabl*   2 
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P<t,a,n) 


rorc«  L«vel,  m 


Force  Level,  n 


figure  1u  -  PLOT  OF  JOINT  PROBABILITIES  P(t,«,n; 
with  data  according  to  table  2 
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k  P(t,in,n) 


force  Level,  m 


Force  Level,  n 


J«l«ur«  11  -  PLOT  OF  JOINT  PROBABILITIES  P(t,«,n; 
with  data  according  to  Table  2 
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P(t,ffl,n) 


Force  Level,    m 


Force  Level f    n 


irigur«   12  -   i'LOT  OF  JOINT   PROBABILITIES   P(t,«,a) 
with  data  according   to  Table  2 
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IV.  EVOLUTION  OF  THE  FORCE  LEVELS 

A.  THE  DETERMINISTIC  MODEL 

The  evolution  of  force  levels  as  a  function  of  time,  X(t)  and  Y(t), 
has  already  been  stated  in  equations  (2.5)  and  (2.6).  The  next  two 
figures  show  the  force  levels  with  different  breakpoints.  It  is  easy  to 
realize  that  introducing  a  nonzero  breakpoint  does  not  change  the  under- 
lying function.  In  other  words,  using  X(t)  as  an  example,  for  both 
figures  the  same  curve  was  used  but  at  the  point  where  the  X-force  reaches 
its  breakpoint,  the  curve  is  "cut."  From  that  point  in  time,  there  are 
no  more  changes  in  the  force  levels.  So  introducting  a  nonzero  break- 
point only  shifts  the  discontinuity  (marked  by  DXl  in  Fig.  13)  up  along 
the  curve  to  DX2. 

The  probability  for  one  side  to  win  is  either  one  or  zero,  because 
it  is  a  deterministic  model.  To  easily  determine  which  side  is  going  to 
win,  a  victory  prediction  condition  can  be  obtained  by  solving  each  force 
level  equation  (2.5)  and  (2.6)  for  the  time  to  reach  its  breakpoint  tx^ 
by  substituting  X(t)=X|^  and  tyj^  by  substituting  Y(t)=y^  .  Then  X  will 
I   win  if  ty.   <  tx.  ,  which  leads  to  the  prediction  condition.  X  will  win 
a  fixed  force  breakpoint  battle  if  and  only  if 

'°''°^   b(x2-x^p)      .         (4.,) 
This  shows  that  given  the  initial  data  one  can  predict  the  outcome 
I   of  the  battle  in  terms  of  force  levels  and  time  until  the  battle 
finishes. 
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^      time  t 


Figure  13  -  FORCE  LEVELS  OVER  TIME 


Deterainistic  Mod«l  with  Breakpoint  x^  =  0 
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^0)Yo 


tirne  t 


Figure  H  -  rORCE  LEVELS  OVER  TIME 


Deterministic  Model  with  Breakpoint  x^  >  0 

Dp' 
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B.  THE  STOCHASTIC  MODEL 

Every  possible  state  (m,n)  in  the  state  space  has  associated  with  it 
a  certain  probability  between  zero  and  one  which  is  a  function  of  time. 
In  order  to  gain  more  insight  into  the  stochastic  process  not  a  single 
realization  of  a  battle  has  to  be  considered,  but  an  average  battle. 
Therefore  the  expected  value  of  the  force  levels  (i.e.  averages)  as  a 
function  of  time  and  the  variances  in  the  force  levels  were  investigated. 

The  straight  forward  way  to  compute  the  expected  force  levels  involves 
the  knowledge  of  the  probability  distribution  P(t,m,n).  Then 

>     "O 
E(M{t)  )=Y^  ^  ^P(t,m,n)  (4.2) 

%p    "bp 

"o    ""o 

E(N(t)  )=y2  "X— ^^^'"^'"^  ^^'^^ 

"bp    %p 

are  the  expected  values  of  the  force  levels  as  a  function  of  time. 

There  are  some  other  ways  to  calculate  the  expected  force  levels, 

one  of  which  will  be  stated  here.  Recall  the  "hybrid-analytical- 

numerical"  method  to  get  the  state  probabilities  (equations  3.5  through 

3.15).  CLARK  [4]  has  also  shown  that  the  i   moment  of,  for  example, 

the  M-force  level  may  be  computed  as 

E(M^t)  )  =  D^'*^  +ZI  ri  D^''^  exp(-(ak+bj)t  )     (4.4) 
^'^   j=l  k=l  ^'^ 

with  for  l^j^mQ  and  Uk^ng 
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0       k 

d(;)  =  5:  .^   cg:S  .    (4.6) 

m=l 

This  again  emphasizes  the  strong  point  of  CLARK' s  solution.  For  a 
given  set  of  battle  parameters  the  coefficients  c'J'J!  have  to  be  computed 
only  once.  Then  with  this  information  and  relatively  small  computa- 
tional effort  not  only  the  state  probabilities  but  also  the  first  and 
second  moment  of  the  force  levels  can  be  computed.  This  determines  the 
variance  in  the  force  levels  at  the  same  time,  e.g. 

Var(M(t)  )  =  E(M^  (t)  )  -  E(M(t))  •  E(M(t))     .    (4.7) 
On  the  other  side,  the  weak  point  is  that  CLARK  [4]  considered  only 
breakpoints  equal  zero. 

Many  authors  have  discussed  one  side's  probability  of  winning  or 
probability  of  winning  conditioned  on  a  certain  number  of  survivors, 
which  always  eliminated  the  parameter  time  by  integrating  from  time  t=0 
to  infinity.  This  may  be  legitimate  to  answer  absolute  (meaning  time 
independent)  questions  about  who  will  win,  but  for  direct  comparisons 
with  the  deterministic  model,  this  author  has  the  feeling  that  the  best 
question  to  ask  regarding  a  winner  is: 

What  is  the  probability  of  one  side  winning  given  the  stochastic  battle 
lasted  as  long  as  the  deterministic  one? 

The  calculation  of  these  probabilities  gives  another  interesting 
probability,  because  given  the  time  t  =  t^  (final  time  of  the  determin- 
istic battle) 
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P(battle  has  not  yet  finished)  = 

1-P(M  wins|t=t^)-P(N  wins|t=t^)         .    (4.8) 

There  has  also  been  work  done  regarding  the  distribution  of  the 
time  to  finish  a  battle.  But  this  is  beyond  the  scope  of  this  thesis 
(SPRINGALL  [9]). 

C.  DIFFERENCES  IN  THE  FORCE  LEVELS 

The  deterministic  model,  especially  the  force  level  equations  (2.5) 

and  (2.6),  describe  a  process  with  continuous  state  parameters  where,  in 

reality,  the  possible  states  are  integers.  Quoting  from  LANCHESTER  [7]: 

Since  the  forces  actually  consist  of  a  finite  number  of  finite  units 
(instead  of  an  infinite  number  of  infinitesimal  units)  the  end  of  the 
curve  must  show  discontinuity,  and  break  off  abruptly  when  the  last 
man  is  reached;  the  law  based  on  averages  evidently  does  not  hold 
rigidly  when  the  numbers  become  small. 

LANCHESTER  suggested  that  his  differential  equations  may  be  good 
approximations  only  as  long  as  the  force  sizes  are  large.  He  also  stated 
that  the  equations  are  based  on  averages,  implying  an  underlying  stochas- 
tic process. 

This  shows  that  there  must  be  a  difference  in  the  force  levels  which 
should  become  significant  when  the  number  of  combatants  is  small.  This 
difference  was  called  bias  by  CLARK  [4]  and  TAYLOR  [10].  It  can  be  shown 
that 

dE(MU))  3  -aE(N(t))+aB^(t)  (4.9) 


and 


^^^^[^^)  =  -bE(M(t))+bB^(t)  (4.10) 


where 
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rriQ  riQ 


and 

"o         '"o 

^m^^^  =  "'bpXl.  P^^'%'"^"'XZ  "iP(t,in,n,^p)  .    (4.12) 
"bp^^  '"bp^l 

The  bias  tenns  B^(t)  and  B  (t)  can  be  interpreted  as  the  expected 
values  of  M(t)  or  N(t)  conditioned  on  the  fact  that  the  battle  has 
already  ended  at  time  t,  for  example, 

B^(t)=E  [M(t))|N(t)=njjp  or  M(t)=mj^pj  (4.13) 

In  other  words,  equation  (4.9)  says  the  expected  casualty  rate  of  the 
M- force  is  proportional  to  the  expected  number  of  survivors  of  the 
N-force  given  neither  of  the  two  forces  has  reached  its  breakpoint. 

Define  the  bias  of  the  X-force  as  Ax(t)  =  E(M(t))-X(t)  and  the  bias 
of  the  Y-force  as  Ay(t)  =  E(N(t))-Y(t) .  Then  using  equations  (2.8)  and 
(2.9)  together  with  (4.9)  and  (4.10)  it  follows  that 

^Ax  =  -aAy+aB^(t)  (4.14) 

^Ay  =  -bAx+bB^(t)  (4.15) 

with  the  initial  conditions  x(0)=)  and  Y(0)=).  This  has  the  solution 

t 
Ax(t)=Vab  y{B^(s)Va7bcosh(Vab' (t-s) 

-Bj^(s)sinh(Vab  (t-s))}ds,    (4.16) 
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and 

t 


Ay(t)=  Vab/  {B|^(s)  Vb/acosh(  Vab(t-s) 


0 

-B^(s)sinh(Vab(t-s))}ds   .         (4.17) 

Since  for  a  fixed,  nonegative  argument  z,  the  cosh(z)  is  always 
greater  than  the  sinh(z),  it  is  easy  to  visualize  that  in  most  of  the 
cases  both  biases  are  positive,  meaning  the  expected  force  levels  of  the 
stochastic  model  are  higher  than  the  deterministic  force  levels.  This 
has  been  shown  by  CLARK  [4]  and  CRAIG  [5]  and  confirmed  by  this  author. 
In  the  rest  of  the  cases  the  winner's  bias  is  negative  or  close  to  zero 
and  the  loser's  bias  is  positive.  Two  examples  are  given  in  Table  3. 
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TABLE  3 

Examples  for  cases  where  one  bias  is  positive  and  the  other 
bias  is  negative. 


A.   Y-force  wins  in  a  fight  to  the  finish 

Xq  =  Yq  =  40  a  =  0.08 


^bp  =  0 

^bp  =  0 

t^   =  15.581 


b  =  0.04 


Ax(t^)  =  3.22 
Ay(t^)  =  -0.23 


B.   Y-force  wins  in  a  fight  with  equal  initial  force  levels, 
but  different  breakpoints. 

x^  =  Yn  =  15  a  =  0.08 

b  =  0.08 


■0  ■  -^0 


bp 

^bp  =  ^ 
t^   =  2.85 


Ax(t^)  =  0.15 
Ay(t^)  =  -0.10 
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An  interesting  point  has  to  be  mentioned  regarding  case  B  of 
Table  3.   Here  the  expected  force  level  of  the  winner  is 
smaller  than  the  expected  force  level  of  the  loser  at  the 
time  a  deterministic   battle  ends. 

To  get  a  better  feeling  for  the  differences  in  the  force 
levels,  the  next  two  figures  show  as  an  example  a  large  spec- 
trum of  force  level  behavior.   The  data  and  notation  is 
described  in  Table  4.   There  are  four  battles  with  four  dif- 
ferent breakpoints  drawn  as  they  evolve  over  time  from  zero 
to  the  time  an  equivalent  deterministic  battle  ends. 
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TABLE    4 

Xq    =    40 

Yq    =    40 

a      =   0.09 

b      =   0.07 

Mi    =    E   M(t) 

Ni    =    E    N(t) 

X^p(i)  =  40-0. 2i 
Y^p(i)  =  40-0. 2i 
for  i  =  1  ,2,3,4 


for  i  =  1,2,3,4 


Where  the  index  i  corresponds  to  the  battle  with 
the  i    breakpoint  force  level: 

X-force  is  always  loser 

Y-force  is  always  winner 
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Figure    1  i?  -    DETERMINISTIC   AND   EXPECTED   FORCE  LEVELS   FOR 
THE  X   -  J?ORCE 
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Flgur©  16  -  DETERMINISTIC  AND  EXPECTED  FORCE  LEVELS  FOR 
THE  Y  -  rORCE 
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CRAIG  [5]  has  formulated  hypotheses  concerning  the  biases  in  the  average 

force  levels  based  on  his  work.  His  hypotheses  were  partially  confirmed, 

but  in  some  cases  they  have  to  be  modified.  Therefore,  another  similar 

set  of  hypotheses  will  be  given  and  supported  by  Fig.  15  and  Fig.  16,  as 

well  as  some  of  the  later  figures. 

H  1)  Given  fixed  initial  force  levels  and  attrition  rate  coefficients, 
as  the  breakpoint  force  levels  increase,  the  numerical  bias  for  the 
loser  decreases.  The  biases  for  the  winner  do  not  show  this  monotone 
trend  except  for  the  case  of  symmetric  parity. 

H  2)  Everything  else  constant,  the  bias  of  the  loser  increases  with 
increasing  initial  force  levels;  this  is  also  true  in  the  symmetric 
parity  case  for  both  forces. 

H  3)  Given  the  initial  force  level  ratio  is  close  to  one  at  the  time 
corresponding  to  the  end  of  the  deterministic  battle,  the  bias  of  the 
loser  is  always  larger  than  the  bias  of  the  winner. 

H  4)  At  the  time  corresponding  to  the  end  of  a  deterministic  battle, 
the  biases  become  larger  as  the  forces  come  closer  to  parity. 

H5)  The  biases  at  times  corresponding  to  less  than  one  half  the 
duration  of  the  deterministic  battle  are  negligible. 

The  case  of  symmetric  parity,  i.e.  equal  initial  force  levels, 

breakpoints  and  equal  attrition  rate  coefficients  seems  to  be  kind  of  a 

"limiting"  case.  For  example,  at  parity  the  biases  of  both  forces  behave 

in  the  same  manner  and  are  equal.  In  Fig.  17  the  biases  at  the  time  a 

deterministic  battle  ends  as  a  function  of  the  initial  force  levels  and 

as  a  function  of  the  breakpoints  are  presented.  It  is  also  another 

verification  for  the  hypotheses  H  1  and  H  2. 
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Figur*  17  -  BIASES  IN  SYMMETRIC  PARITY 
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Fi«ur«  18    -   BIASES  WITH   DIl-TERENT    INITIAL  FORCE  LEVELS 
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Fi«ur«  19  -  BIASES  WITH  DIFFERENT  ATTRITION  RATE 
COEFFICIENTS 
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Considering  the  changes  in  magnitude  of  the  biases  for  battles  like 
in  Fig.  15,  16  and  17,  CRAIG  [51  came  to  the  conclusion  that  when  the 
forces  are  closer  to  parity,  the  biases  at  the  deterministic  battle's 
end  increase.  Several  battles  were  fought  starting  with  symmetric  parity 
and  then  varying  the  force  levels  and  the  attrition  rate  coefficients. 
Sample  results  are  shown  in  Fig.  18  and  Fig.  19.  In  Fig.  18,  the  initial 
force  levels  were  changed  giving  the  Y-force  ten  and  twenty  percent 
higher  initial  force  levels.  The  biases  at  the  end  of  the  deterministic 
battle  are  plotted  as  a  function  of  the  breakpoint  force  level  ratio 

^  ^   ^bp/^0- 

Fig.  18  supports  CRAIG' s  hypothesis,  as  does  Fig.  19.  Here  not  the 
initial  force  levels  but  the  attrition  rate  coefficients  were  changed  in 
order  to  deviate  from  symmetric  parity.  The  last  way  to  deviate  from 
symmetric  parity  is  a  case  where  the  deterministic  model  gives  equal 
answers  for  different  battles.  The  data  and  the  results  are  shown  in 
Table  5. 
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TABLE  5 
Equal  initial  force  level  battle  with  non  equal   breakpoints 


Xq  =  15        Yq  =  15        a  =  b  =  f  =  0.08 


X^p  Y^^p  Ax(t^)  Ay(t^) 

12  12  1.08  1.08 

12  9  0.71  0.63 

12  6  0.15  -0.10 

12  3  0.68  0.58 

9  9  1.62  1.62 

9  6  1.14  0.93 

9  3  1.06  0.75 

6  6  2.15  2.15 

6  3  1.67  1.22 

3  3  2.84  2.84 
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The  X-force  level  at  time  t^  is  the  same  as  the  Y-force  level.  Thus, 
the  expected  force  level  for  the  winner  (i.e.  the  Y-force)  is  smaller  than 
the  expected  force  level  of  the  loser  (the  X-force),  because  Ax(t^)  is 
always  larger  than  Ay(tr).  This  is  easy  to  see  when  one  remembers  the 
way  the  expected  force  level  is  computed  (equation  (4.3)).  Since  the 
Y-force  has  a  lower  breakpoint,  there  are  states  (m,n)  possible  where 
n.  <n<m.  .  Apparently  these  states  have  a  nonzero  probability  associated 
with  them,  which  decreases  the  expected  force  level  below  the  expected 
value  for  the  X-force.  This  might  be  a  starting  point  for  further 
studies. 

D.  VARIABILITY  IN  THE  FORCE  LEVELS 

Naturally  in  the  deterministic  case  there  does  not  exist  any  variabil- 
ity in  the  force  levels.  On  the  other  hand,  for  the  stochastic  model,  the 
variance  in  the  force  levels  as  a  function  of  time  is  a  measure  of  the 
dispersion  of  the  number  of  survivors  about  their  mean  value. 

CLARK  [4]  has  hypothesized  two  different  types  of  behavior  for  the 
variance  in  the  force  levels,  which  are  shown  in  Fig.  20  for  the  data 
presented  in  Table  2.  The  first  type  of  behavior  is  that  of  the  variance 
for  the  N-force,  i.e.  the  variance  increases  monotonely  as  a  function  of 
time  and  is  asymptotic  to  a  limiting  value.  It  was  found  that  this  type 
of  behavior  occurs  when  the  side  is  going  to  win  and  for  the  case  of 
symmetric  parity.  The  second  type  of  behavior  shown  is  the  variance  of 
the  M-force,  Var(M(t)),  as  a  function  of  time.  This  increases  to  a 
maximum  value  then  decreases  asymptotically  to  a  limiting  value.  This 
type  of  behavior  is  associated  with  the  loser  of  the  battle. 


63 


The  variance  in  the  force  levels  is  a  function  of  the  initial  and 

breakpoint  force  levels,  the  battle  time  and  the  attrition  rate 

coefficients.  Unfortunately  one  does  not  know  what  this  dependence  is. 

Based  on  many  numerical  results  a  set  of  hypotheses  will  be  stated  and 

the  next  figures  will  support  them. 

H  1)  Given  fixed  initial  force  levels  and  attrition  rate  coefficients, 
as  the  breakpoint  force  levels  increase  the  variance  of  the  force  levels 
decrease. 

H  2)  Everything  else  constant  the  variance  of  the  loser's  force  level 
increases  with  increasing  initial  force  levels.  This  is  also  true  for 
both  variances  in  the  case  of  symmetric  parity. 

H  3)  Given  the  initial  force  ratio  is  close  to  one  at  the  time 

corresponding  to  the  end  of  the  deterministic  battle  the  variance 

of  the  loser's  force  level  is  always  smaller  than  the  winner's  variance. 

H  4)  At  the  time  corresponding  to  the  end  of  the  deterministic  battle 
the  variance  in  the  loser's  force  level  increases  as  the  forces  come 
closer  to  parity.  This  trend  is  not  true  for  the  variance  of  the  win- 
ner, except  for  the  case  of  symmetric  parity. 

It  was  the  intention  of  this  author  that  the  set  of  cases  to  illus- 
trate the  hypotheses  are  the  same  as  in  the  illustrations  (Fig.  17,  18 
and  19)  of  the  hypotheses  about  biases.  So,  Fig.  21  shows  the  force 
level  variances  for  the  different  initial  force  levels  and  different 
breakpoints  for  the  case  of  symmetric  parity.  In  Fig.  22  the  variances 
for  a  battle  with  equal  attrition  rate  coefficients  but  varying  initial 
force  levels  show  that  H  4  is  only  true  for  the  loser.  This  point  is 
emphasized  by  Fig.  23,  where  with  constant  and  equal  initial  force  levels 
the  attrition  rate  coefficients  were  varied. 
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rifur*  20  -  FORCE  LEVEL  VARIANCES  OVER  TIME 
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Jrisup«  21    -   VARIANCES    IN  SYMMETRIC   PARITY 
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Figure  22  -  VARIANCES  WITH  DIFFERENT  INITIAL  FORCE  LEVELS 
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V.  COMPUTATIONAL  ASPECTS 

A.  USING  THE  NUMERICAL  SOLUTION 

For  most  of  the  numerical  work,  the  state  probabilities  have  been 
obtained  using  the  fourth-order  RUNGE-KUTTA  method.  The  accuracy  of  the 
results  was  increased  by  substituting  the  available  analytical  results 
for  the  Regions  I,  II  and  III. 

The  disadvantage  is  that  this  solution  method  needs  a  lot  of  CPU-time. 
For  the  battle  in  Example  1,  with  the  data  given  in  Table  1,  the  calcula- 
tion of  the  state  probabilities,  expected  force  levels  and  variances  led 
to  a  CPU-time  on  an  IBM-360  computer  of  almost  90  minutes  with  a  time 
step  size  of  0.05  minutes.  That  was  the  reason  why,  in  the  examples  the 
attrition  rate  coefficients  were  increased  by  one  magnitude,  which  brought 
the  computer  usage  down  to  a  CPU-time  of  around  12  minutes  with  the  same 
time  step  size  of  0.05  minutes.  The  sum  of  the  state  probabilities  at 
every  time  step  was  used  as  a  measure  of  accuracy.  Surprisingly,  its 
deviation  from  1.0  was  always  in  the  fourth  or  fifth  decimal,  which  proves 
the  robustness  of  the  RUNGE-KUTTA  method.  Also,  the  fact  that  for  the 
Regions  I,  II  and  III  analytical  solutions  were  substituted  at  each  time 
step  did  not  change  the  final  outcome  considerably.  It  was  found  that 
without  the  analytical  partial  solutions  the  sum  of  probabilities  tended 
to  be  slightly  higher,  changing  the  sum  of  probabilities  from,  for 
example,  0.99998  to  1.00003. 
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B.  THE  EQUAL  ATTRITION  RATE  COEFFICIENT  SOLUTION 

After  the  development  of  the  equal  attrition  rate  coefficient  solution 
(EARCS)  outlined  in  Ch.  III.C,  it  was  implemented  on  an  IBM-360  comput- 
er. The  CPU-time  for  the  calculation  of  the  state  probabilities,  expected 
force  levels  and  variances  for  a  given  point  in  time  always  stayed  below 
20  seconds,  which  emphasizes  its  computational  advantage. 

Further  research  showed  two  shortcomings  of  the  EARCS,  which  are 
easily  overlooked.  The  coefficients  C(m,n)  vary  over  a  wide  range  start- 
ing at  1.0  and  depending  on  the  initial  force  levels. 
For  example  for  a  battle  with 

mQ  =  20  mj^p  =  15 

"O  =  ^0  "bp  =  20 

f  =  0.08        t^  =  1.60995 

C(20,40)  =  1.0    but   C(16,21)  =  4.9297077276x10"^^. 
Also,  the  binomial  coefficients  in  equations  (3.23)  and  3.24)  vary  over  a 
wide  range  starting  at  1.0  to  1352078  for  the  above  example.  This  shows 
that  the  capacity  the  computer,  in  terms  of  the  number  of  significant 
digits,  is  able  to  carry  limits  the  implementation  of  this  solution  in 
the  present  form. 

Another  reason  why  this  solution  is  not  the  end  of  the  numerical 
problems  is  the  fact  that  in  the  given  form,  equations  (3.23)  and  (3.24) 
require  computation  of  an  alternate  sum  consisting  of  terms  whose  factors 
are  relatively  large  binomial  coefficients  and  small  numbers  between  zero 
and  one.  This  produces  truncation  errors,  which  yield  nonsensical  results 
like  negative  variances  and  sums  of  probabilities  greater  than  one. 
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VI.  CONCLUSIONS 

Throughout  this  thesis,  only  LANCHESTER-type  square-law  attrition 
with  fixed  initial  force  levels  and  fixed  attrition  rate  coefficients 
has  been  considered  for  a  deterministic  and  a  stochastic  version.  The 
stochastic  version  required  much  more  computational  effort.  So,  given 
the  need  for  an  analytical  model  as  opposed  to  the  use  of  simulation, 
there  is  not  much  to  gain  from  the  application  of  a  stochastic  model 
when  the  force  levels  are  large  and  the  forces  are  not  near  parity.  In 
these  cases  the  deterministic  version  essentially  produces  the  same 
results,  at  least  in  qualitative  terms.  For  smaller  force  levels  or 
forces  near  parity,  the  stochastic  version  may  be  helpful  to  get  more 
information  about  the  dynamics  of  combat. 

Also,  there  exist  one  case  where  the  deterministic  version  cannot 
differentiate  between  several  different  battles.  That  is  the  case  of 
equal  initial  forces,  equal  attrition  rate  coefficients,  but  different 
breakpoint  force  levels.  This  also  rectifies  the  further  development  of 
the  equal-attrition-rate-coefficient-solution  (EARCS). 

Using  an  analytical  model,  the  discussed  way  of  introducing  random- 
ness into  the  model  seems  not  to  be  very  enlightening.  Therefore,  it  is 
suggested  that  another  way  to  include  random  effects  should  be  explored. 
The  author's  opinion  is  that  working  with  attrition  rate  coefficients 
which  are  random  variables  seems  more  promissing  to  gain  insight  into 
the  dynamics  of  combat.  Further  down  the  line  there  should  be  some 
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consideration  on  the  usage  of  combinations  of  the  possibilities  to 
include  random  effects  as  given  in  Chapter  II. B. 
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APPENDIX  A 


1 
2 
3 

« 
S 

6 

7 

8 

9 

10 

11 

12 

13 

lU 

15 

16 

17 

16 

19 

20 

21 

22 

23 

2M 

25 

26 

27 

28 

29 

30 

31 

32 

33 

311 

35 

36 

37 

36 

39 

UO 

til 

42 

Hi 

W 

MS 

V6 

«7 

48 

49 

SO 


THIS  PROGRAM  COMPUTES  THE  STATE  PROBABILITIES  USING  THE  FOURTH 
ORDER  RUNGE-KUTTR  METHOD.  PLOTS  OF  THE  EXPECTED  FORCE  LEVELS  AND 
VARIANCES  AS  HELL  AS   3-D-PLOTS  FOR  THE  STATE  PROBABILITIES  ARE 
OBTAINED. 

DIMENSION  PI  (33,33),  P2(33,33).  T3D  (20) 

DIMENSION  XMT  (500)  ,  XNT  (500)  ,  OT  (500)  ,  VM  (500)  ,  VN  (500)  , 
1  VMB(551),  VNB(551),  VQ(551),  VARXtSSl),  VARY  (551),  TIM(551) 

DIMENSION  EST  (313)  .0(313)  .F  (2)  .SIZE  (2)  .  KX  (100)  ,  KT  (100)  .  WK  (41 ,  41 ,  3) 
1        ,  PC  (m  ,11)  ,CST  (mi  .CETERM  (41) 

LOGICALnI  I0N(41,41) 

DIMENSION  DETERM(313) 

REflLxS  TTL  (12) /12X'         V 

REAL  Kl,  K2,  K3,  KM 

CALL  ERRSET (208.600.-1.1) 

F(1)«0. 

F(2)«=0. 

LlNES-0 

SIZE  (l)-6. 

SIZE (2) -6. 

OFLAG  -  0. 

NKXT>100 

READ  (5,211 

READ  (5,22) 

READ  (5.22) 

READ  (5.21) 

READ  (5,22) 

TBDd) 

IF  (N3D 


MBP,  MO,  NBP.  NO 
AA.  BB 

H,  FTIM 
N3D 

EPSQ,  EPSPTT 
9993. 
.LT.  1) 


READ  (5, 22) 
25  CONTINUE 


GO  TO  25 
(T3D(I),  1  =  1. N3D) 


N3D-NUMBER  OF  3D-PL0TS,   H=SIZE  OF  TIME  STEP 
FTIHE-FINALTIME.  EPSQ. EPSPTT   ARE  ZERO  LEVELS 
T3D(I)-  FRACTION  OF  FTIME  WHEN  TO  PLOT  3D 

00  10  I-1.N30 
T30(I)-FTIM«T3D(I) 
10  CONTINUE 

21  FORMAT  (1615) 

22  FORMAT  (8F10. 5) 

WRITE  (6,28)  MBP.  MO.  NBP.  NO 

28  FORMAT  (//  5X,  'M  S  «  N  S',  5X.  416  /) 
HRITE(6.29)  AA,  BB 

29  FORMAT  (/  5X,  'A  I    B',  9X.  2F10.3  /) 
MR1TE(6,30)   H,  FTlM 

30  FORMAT!/  5X,  'INITIAL  H  *  FINAL  TIME  (•  LOOPS  x  H) '.  5X,  2F10.3/) 
WRITE  (6.31)  N3D 

31  FORMAT  (//  5X,  '«  OF  3-D  PLOTS',  2X.  16  /) 
IF  (N3D  .IT.  1)  GO  TO  36 
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51  WRITE  (6,214)  EPSQ.EPSPTT 

52  2q  FORMAT  (/  5X. 'EPSO= * . F 10. 5. 5X.  •EPSPTT= ' , FIO. 5) 

53  WRITE  (6, 32)  (TSOd).  1  =  1. N3D) 

51  32  FORMfiT(/  5X.  'AT  TIME'.  2X.  IOFIO.3  /) 

55  36  CONTINUE 

56  MOl  »  MO  +  1 

57  NOl  »  NO  ♦  1 

58  NS  ■  NOl  -  NBP 

59  MS  -  MOl  -  MBP 

60  MSL  -  MS  ♦  1 

61  NSL  «  NS   +1 

62  BMO  >  BB  M  MO 

63  ANO  -  RA  M  NO 

64  QANB  >  ANO  /  BB 

65  OBMA  '   BMO  /  AA 

66  BAT  =  -  (BMO  *  ANO) 

67  DO  47  Ml  =  I.MS 

68  VM  (Ml)  '    MBP  ♦  Ml  -  1 

69  47  CONTINUE 

70  DO  48  Nl  =  l.NS 

71  VN(Nn  =  NBP  ♦  Nl  -  1 

72  48  CONTINUE 

73  L  -  0 

74  N3  -  1 

75  TPL0T  -  T30(N3) 

76  TIME  =.  0. 

77  50  CONTINUE 

78  IF  (TIME  .GT.  FTIM)  GO  TO  210 

79  52  CONTINUE 

80  TIME  -  H  X  L 

81  L  -  L  ■►  1 

82  LB  »  0 

83  TIM(L)  =  TIME 

84  IF  (  L  .EO.  2)  GO  TO  53 

85  IF  (TIME  .LT.  TPLOT)  GO  TO  60 

86  C 

87  C  SET  NEXT  PLOT  TIME 

88  C 

89  N3  »  N3  ♦  1 

90  TPLOT  »  T3D  (N3) 

91  IF  (N3.GT.N3D)  TPL0T  =  9999. 

92  53  CONTINUE 

93  C 

94  C  LB  =  L  WHEN  TIME  TO  PLOT 

95  C 

96  LB  -  L 

97  EBT  -  EXP(BB  «  TIME)  -  1. 

98  EAT  -  EXP(AA  »«  TIME)  -  1. 

99  AE  '   QANB  x  EBT 

100  BE  "    OBMA  «  EAT 


74 


101  BflE  «  EXP(BflT  w  TIME) 

102  FJ  -  1. 

103  JJ  -  1 
1011  C 

105  C  NO  VALUE  OF  J  =  ZERO  fiT  MO. NO 

106  C 

107  FK  -  1. 

108  KK  -  1 

109  QFK-1. 

110  QFJ-1. 
HI  C 

112  C  NO  VALUE  OF  K  =  ZERO  AT  MO, NO 

113  C 

im  60  CONTINUE 

lis  IF  (L  .EQ.  n  GO  TO  500 

116  00=^0. 

117  DO  400  Nl  =  l.NS 

118  N  «  NOl  -  Nl 

119  NL  -  NSL  -  Nl 

120  NLl  -  NL  -»  1 

121  ANJ  «  AA  V  N 

122  C 

123  C 

1211  00  300  Ml  »  l.HS 

125  H  «  MOl  -  Ml 

126  ML  *  MSL  -  Ml 

127  ML  1  -  ML  ♦  1 

128  BHI  >  BB  M  M 

129  ABA  -  ANJ  *   BMI 

130  IF  (M  .EQ.  MO)  GO  TO  502 

131  IF  (M  .EQ.  MBP)  GO  TO  504 

132  IF  (N  .EQ.  NO)  GO  TO  506 

133  IF  (N  .EQ.  NBP)  GO  TO  507 
1311  C 

135  C  DEFAULTS  TO  ALL  INSIDE  POINTS 

136  C 

137  Kl  »  ANJwP2  (MLl.NL)  *  BMl«P2  (ML. NL 1 )  -  RBA«P2  (ML.NL) 

138  PT  -  ANJmO.Sx  (PI  (MLl.NL) +P2  (MLl.ND)  ♦  BMI«0.5  »<  (PI  (ML. NLl)  ♦ 

139  1  P2(ML.NLn) 

140  K2   -    PT    -    ABA   w     (P2(ML,NL)     ♦    HmO.SxKD 
im  K3   ■   PT    -    ABA   H    (P2(ML.NL)     *    H»<0.5»<K2) 

142  K4    ■   flNJxPl  (MLl.NL)     ♦    BMI«P1  (ML.  NL  1 )     -    RBRm  (P2  (ML.  NL)     +    H»«K3) 

ll|3  PI  (ML.NL)     -    P2(ML.NL)    +     (H/6.0)     »<     {K1*2.0mK2    ♦    2.0mK3    ♦    Ki4) 

mV  GO   TO    200 

145  502    IF    (N    .EQ.    NO)    GO   TO   503 

146  IF    (N    .EQ.    NBP)    GO   TO   507 

147  C 

148  C  M=MO,    N-=N0,NBP 

149  C 

150  Kl    ■   BMI»«P2(ML,NLn     -    RBR«P2  (ML.  NL) 
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151 

PT  -  BHixO.S  «  (PMML.NLl)  ♦  P2(ML.NLn) 

152 

K2  -  PT  -  flBfi  X     (P2(ML.NL)  *  HkKImO.5) 

153 

K3  -  PT  -  flBfl  H  (P2(ML,NL)  ♦  HvKS^O.Sl 

154 

Ky  =•  BMl  X  PMML.NLl)  -  flBfl  «  (P2(ML.NL)  + 

155 

PICML.NL)  =•  P2(ML.NL)  *  (H/6.)  «  (K1-^2.0hK 

156 

IF  (LB  .EQ.  L)  GO  TO  601 

157 

GO  TO  200 

158 

C 

159 

C 

M'-MO,  N=NO 

160 

C 

161 

503 

CONTINUE 

162 

P1(ML,NL)  =  EXP  (-flBfl  H  TIME) 

163 

HRITE(6,102)  L.  N,  ML.NL,  PMML.NL) 

164 

GO  TO  200 

les 

504 

IF  (N  .EG.  NBP)  GO  TO  505 

166 

C 

167 

C 

168 

C 

M=MBP,  N-=  NBP 

169 

C 

170 

Kl  ■  flNJ  w  P2(ML1.NL1 

171 

K2  -  fiNJ  n    0.5  •«  (PMMLl.NL)  ♦  P2(ML1,NL)) 

172 

K3  -  K2 

173 

K4  -  flNJ  H  PI  (MLl.NL) 

174 

PI  (ML.NL)  =«  P2  (ML.NL)  ♦  (H/6.0)  x  (Kl  +  2. 

175 

GO  TO  200 

176 

c 

177 

c 

M=MBP.  N=NBP 

178 

c 

179 

505 

CONTINUE 

180 

PI  (ML.NL)  -0.0 

181 

GO  TO  200 

182 

c 

183 

c 

N»NO.  M-=M8P,  MO 

184 

c 

185 

506 

Kl  «  AN  JXP2  (MLl.NL)  -  flBfl»«P2  (ML.  NL) 

186 

PT  -  flNJ  «  0.5  X  (PI  (MLl.NL)  ♦  P2(ML1.NL)) 

187 

K2  -  PT  -  flBfl  «  (P2(ML.NL)  +  H«0.5»<Kn 

188 

K3  -  PT  -  flBfl  ♦•  (P2  (ML.NL)  ♦  H«0.5»«K2) 

189 

K4  •  flNJ  M   PI  (MLl.NL)  -  flBfl  n     (P2  (ML.NL)  ■► 

190 

PI  (ML.NL)  »  P2  (ML.NL)  +  (H/6.0)  h     (Kl  >  2. 

191 

IF  (LB  .EQ.  L)  GO  TO  602 

192 

GO  TO  200 

193 

c 

194 

c 

N=NBP,  t,    N=NBP.  M=MO 

195 

c 

196 

507 

CONTINUE 

197 

Kl  ■  BMI  n    P2(ML.NL1) 

198 

K2  -  BMI  K  0.5  H  (Pl(ML.NLl)  ♦  P2(ML.NL1)) 

199 

K3  »  K2 

200 

K4  -  BMI  m    PI (ML.NLl) 

HxK3) 
?  ♦  2.0xK3  ♦  K4) 


2.0«K2  ■»  2.0xK3  ♦  Kil) 


HxKS) 
0hK2  *  2.0wK3  ♦  K4) 
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201  PI  (HL.NL)  =  P2(ML.NH  +  (H/S.O)  «  (Kl  *  2.0»«K2  *  ?.0»<K3  ♦  KM) 

202  GO  TO  200 

203  C 

20«1  C  ALSO  M=M0.  N-=NBP.N0 

205  C 

206  601  CONTINUE 

207  IF  ((BE.LT. 1.) .RND.  (KK.GT.9) )  GO  TO  255 

208  C 

209  C     COMPUTE  PTT  IF  BE  >  0  4  KK  <  10 

210  C 

211  PTT=BflE 

212  DO  250  IN0=1.KK 

213  XIN0>IN0 

2111  PTT-PTTkBE/XIND 

215  250  CONTINUE 

216  GO  TO  260 

217  255  PTT-0. 

218  260  CONTINUE 

219  KK  ■  KK  ♦  1 

220  lF((fiBS(Pl  (ML.NL)  -PTT)) .LE.EPSPTT)  GOTO  199 

221  GO  TO  700 

222  C 

223  C  M^MBP,  N=«N0 
2211  C 

225  602  CONTINUE 

226  IF C  (fiE.LT.l.) .AND.  (JJ.GT.9) )  GO  TO  270 

227  C 

228  C 

229  C     COMPUTE  PTT  IF  flE  >  0  4  JJ  <  10 

230  C 

231  PTT-BflE 

232  DO  265  INO^^I^JJ 

233  XIND>INO 

2311  PTT-PTT-fiE/XIND 

235  265  CONTINUE 

236  GO  TO  275 

237  270  PTT=0. 

238  275  CONTINUE 

239  JJ  -  JJ  ♦  1 

2U0  IFKfiBSfPl  (ML.NL)    -PTT)  )  .LE.  EPSPTT)    GO   TO    199 

2m  700   CONTINUE 

2(12  C 

2H3  C     TIME  REDUCED  1/2 

244  C 

245  L  •  L  -  1 

246  HRITE  (6,701)  TIME,  M,  N.  H 

247  701  FORMAT  (//  5X,  'H  VALUE  IS  REDUCED  BY  HALF  AT  :  TIME  =  '.  F8.3. 

248  1   •   M  »  ',  15,  '   N  =  '.  15,  '   FROM  '.  F8.3  //) 

249  H  -  0.5  M  H 

250  WRITE  (6, 102)  L.M.N.KK,  PTT.  Pi  (ML. NL) , QFK.  QFJ.QO.QFLAG 
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251  102  FORMAT  (/  2X,  415,  7F12.5) 

252  QFLflG=0. 

253  GO  TO  52 

25<1       199  PTT  =  P1  (ML.NL) 

255  200  CONTINUE 

256  QQ=QQ*P1  IML.NL) 

257  300  CONTINUE 

258  >400  CONTINUE 

259  QQl=flBS(QO-l.) 

260  IF(QQl.LE.EPSQ)  GO  TO  401 

261  QFLRG-1. 

262  GO  TO  700 

263  101  CONTINUE 

264  C 

265  C      COMPUTE  VflRX  IT)  ,  VARY  (T)  .  N3P1R  (T)  .MBflR  (T)  ,Q  (T) 

266  C 

267  SMB  =  0. 

268  SNB  =  0. 

269  SO  =  0. 

270  SM2  '   0. 

271  SN2  -  0. 

272  DO  415  Ml  =  I.MS 

273  RM  -  VM  (Ml) 
27i|  RM2  =  RMhRM 

275  DO  410  Nl  =  l.NS 

276  RN  "    VN  (Nil 

277  RN2  =  RN«RN 

278  PT  =  PI  (Ml.Nl) 

279  P2(M1,N1)  =  PT 

280  SMB  =  SMB  +  RM«PT 

281  SNB  =«  SNB  ■►  RN»<PT 

282  SO  =  SO  ♦  PT 

283  SM2  =  SM2  ♦  RM2  x  PT 
28l|  SN2  «  SN2  ♦  RN2  «  PT 

285  410  CONTINUE 

286  415  CONTINUE 

287  VMB  (L)  =  SMB 

288  VNB  (L)  ^    SNB 

289  VQ(L)  =>  SQ 

290  VflRX  (L)  =  SM2  -  SMB»<5MB 

291  VARY  (L)  =  SN2  -  SNB»«SNB 

292  C 

293  C     RETURN  TO  MAIN  LOOP  (50)  IF  NOT  TIME  TO  PRINT 
2911  C 

295  IF  (N30  .LT.  1)  GO  TO  430 

296  IF  (LB.NE.L)  GO  TO  430 

297  IF(LB.EQ.2)  GO  TO  430 

298  C 

299  C     ADJUST  X.Y.  VECTORS  FOR  PLOT 

300  C 
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301 

DO  15  K-1,MS 

302 

15 

EST  (K) =K-1 

303 

DO  16  K=1,NS 

304 

DETERM(K)=K-1 

305 

16 

CONTINUE 

306 

DO  115  1-1.41 

307 

CST(n-I 

308 

CETERMH)-! 

309 

DO  115  J* 1.41 

310 

PC(l,J)-0.0 

311 

115 

CONTINUE 

312 

MP«41 

313 

NP-41 

3m 

flLP"15. 

315 

BETfi=30. 

316 

PMflX  =  Pl  (1. 1) 

317 

DO  116  1=1, NS 

318 

HR1TE(6.199)  I. PI  (1  .1)  .PI  (1.1  ) 

319 

199 

FORMfiTC  ',I5.10X.2F15.5) 

320 

116 

CONTINUE 

321 

DO  18  1=1. MS 

322 

DO  18  J=1.NS 

323 

IFIPI  (1,J)  .GT.PMftX)  PMflX  =  Pl(I,J) 

32«1 

18 

CONTINUE 

325 

C0NST=8.0/PHflX 

326 

00  19  1=1. MS 

327 

N=I+MBP 

328 

00  19  J»1.NS 

329 

M=J*NBP 

330 

PI  (I.J)=C0NSTkP1  (I.J) 

331 

PC(N.M)=P1  (l.J) 

332 

19 

CONTINUE 

333 

CALL  PLT301  (CST. MP. CETERM. NP. PC. ALP. BETA. F. TTL. SIZE. WK, IDN. KX.  KT.N 

334 

IKXT. LINES) 

335    C 

336 

430 

CONTINUE 

337    ( 

338    ( 

339    ( 

PLOT  HERE 

3U0    ( 

RETURN  TO  MAIN  LOOP 

341    ( 

3112 

GO  TO  50 

3143 

500 

CONTINUE 

31414 

DO  501  I  =  I.MS 

345 

DO  501  J  =  l.NS 

346 

PI  (I.J)  =  0.0 

347 

501 

CONTINUE 

348 

PI  (MS.NS)  =  1.0 

349 

GO  TO  401 

350    ( 
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351 

C 

FINAL  TIME  REACHED 

352 

C 

353 

210 

CONTINUE 

35^ 

WRITE  (6.211)  TIME,  FTIM,  L 

355 

211 

FORMAT  (//  5X.  'COMPUTED  TIME  '.  Fl 

356 

1  F10.3.  5X.  ••  OF  LOOPS  TO  REACH  F 

357 

DO  215  I  -l.L 

358 

WRITE(6,213)  VMBII),  VNB(I),  VQtl) 

359 

213 

F0RMAT(2X,  12F10.5) 

360 

215 

CONTINUE 

361 

650 

FORMAT  CI        •) 

362 

WRITE  (6,650) 

363 

CALL  PLOTP"(TIM,VMB.L,0) 

3611 

WRITE  (6,650) 

365 

CALL  PLOTP  (TIM, VRRX.L.O) 

366 

WRITE  (6,650) 

367 

CALL  PLOTP  (TIM.VNB,L,0) 

368 

WRITE  (6.650) 

369 

CALL  PLOTP  (TIM, VART.L.O) 

370 

STOP 

371 

DEBUG  SUBCHK 

372 

END 

F10.3.  5X. 

'INPUT  FINAL  TIME 

FINAL  TIME 

'.  15  //) 

I)  ,  VARX  (I)  . 

,  VARY  (I)  ,  TIM(I) 
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APPENDIX  B 


I 

2 

9 
4 
5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 
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THIS  PROGRAM  CflLCULflTES  THE  STATE  PROBABILITIES  FOR  THE  EQUAL 

ATTRITION  RATE  COEFFICIENT  SOLUTION  (EfiRCS) . 

fl   3-D  PLOT  IS  PRODUCED  USING  THE  VERSATEC  PLOTTER. 

IMPLICIT  REALkB  (A-H.O-Z) 

CALL  ERRSET  (208.256.10. 1) 
DIMENSION  C  (50.50)  .M  (50)  .N  (50)  .F  (50)  .PTMN  (50.50) 


m.3)  .X  (WD  .T  (Ml)  ,p  (m.m) 


V 


REflLxM  SIZE  (2)  .FL  (2)  .HK  (Ml 
DIMENSION  KX  (100) ,KT (100) 

logicalmi  ]DN(m.m) 

REALn8  TTL(12)/12«' 
READ  (5, lOOlHO.MBP.NO.NBP 
READ  (5, 101) A 
READ (5, 102) TIME 
READ  (5,103)  (F  (I)  .  1  =  1.41) 


F(I)-1/1-FACT0RIAL.   DONE  TO  SPEED  UP  THE  PROGRAM 


WRITE (6,B02)M0.MBP.N0,NBP.A,TIME 
802  FORMAT  ('   *, 415. 2F 10. 5) 
MD=MO-MBP 
ND=NO-NBP 
M01=M0-^1 
N01=N0*1 
MBP1=MBP*1 
NBP1=NBP*1 
MBP2=MBP*2 
NBP2=NBP*2 
RMO=Ma 
RNO^NO 
DO  10  I»l,41 
DO  10  Jsl,41 
C(l,  J)-0.0 
PTMNd,  J)=0.0 
P(I.J)=0.0 
10   CONTINUE 

C  (MOl.NOn  =1.0 

DO  20  1=2, MO 

J=MBP-»I 

C(J,N01)=  RNO»<«  (MD-I*n 

20  CONTINUE 

DO  21  1=2, NO 

J=NBP*I 

C(M01,J)=  RMOxK  (ND-nn 

21  CONTINUE 

DO  22  1=2, MO 
MC=M01-I 
MM=M01-I*1 
MPLUS=MM*1 
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51 

DO  22  J=2,ND 

52 

NC=N01-J 

53 

NN=N01-J*1 

54 

NPLUS=NN*1 

55 

C  (MM.NN)=NC»«   C  (HPLUS.NN)  ♦HC*.   CCMM.NPLUS) 

56 

22 

CONTINUE 

57 

C 

58 

C 

END 

OF  COEFFICIENT  CflLCULflTION 

59 

C 

60 

n--fl«TIME 

61 

DO  25  I=MBP2.M01 

62 

DO  25  J=NBP2,N01 

63 

IJ=1*J 

611 

K=M0*N0-IJ  *2 

65 

F2"  (1.0-DEXP  (Fl)  )  *>«K 

66 

KK=K*1 

67 

lJ=lJ-2 

68 

F3=DEXP(F1»«IJ) 

69 

PTMNd,  J)=F2«F3«C  (I,  J)  ».F  (KK) 

70 

25 

CONTINUE 

71 

100 

FORMAT  (415) 

72 

101 

FORMAT  (F10.5) 

73 

102 

FORMAT  (F10.5) 

7i| 

103 

FORMAT  (D17.1 1) 

75 

C 

76 

C 

FOR 

NBP<N<NO 

77 

C 

78 

DO  30  11=1. NO 

79 

NFORCE  =  NBP-^II 

80 

J=M0*N0-MBP-1-NF0RCE 

81 

FACT=MBP*1.0+NF0RCE 

82 

SUM1=  (1.0-DEXP  (F  IMPACT))  /  (  n»<FACT) 

83 

SUHINT=SUM1 

84 

DO  301  K-l.J 

85 

FACT»FACT*1.0 

86 

flDDFAC= (-1.0)  ««K 

87 

FRACTN=  (1.0-DEXP (FIkFACT) )/(  A«FflCT) 

88 

C0MBT=1.0 

89 

DO  302  KJ=1,K 

90 

RKJ«KJ 

91 

C0MBI=(J-RKJ*1.0)/RKJ 

92 

COMBT-COMBT«COMBI 

93 

302 

CONTINUE 

911 

SUMlNT=SUMINT*ADOFflC«FRflCTNKCOMBT 

95 

301 

CONTINUE 

96 

NMBP=J*1 

97 

NBOUND=NFORCE*I 

98 

PFAC=AxNFORCE«F (NMBP) kC (MBP2. NBOUND) 

99 

PTHN(MBPl.NBOUND)  =PFRCmSUM1NT 

100 

30 

CONTINUE 
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101 

C 

102 

C 

FOR 

MBP<M<Ma 

103 

c 

lOU 

DO  ^0  11=1. HD 

105 

MF0RCE=MBP+11 

106 

J»M0*N0-NBP-1-MF0RCE 

107 

FflCT  =  NBP-^1.0*HFORCE 

108 

SUM1=  (1.0-DEXP(Fl«FflCT)  )/(  fl^FflCT) 

109 

SUMINT=SUM1 

110 

DO  401  K=1.J 

111 

FflCT=FfiCT*1.0 

112 

flDDFflC=  (-1.0) ««K 

113 

FRflCTN=  (1.0-DEXP(Fl«FfiCT))/  (  fl^FRCT) 

114 

COMBT=il.O 

115 

DO  402  KJ=1,K 

116 

RKJ=KJ 

117 

COMBI=  (J-RKJ+1.01 /RKJ 

118 

C0MBT»C0MBThC0MB1 

119 

402 

CONTINUE 

120 

SUMlNT=SUMlNT+RDDFflCxFRflCTNMCOMBT 

121 

401 

CONTINUE 

122 

MNBP=J-^1 

123 

MBOUND^MFORCE*! 

124 

PFflC  =  R»«MFORCE»<F  (MNBP)  »<C  (MB0UND.NBP2) 

125 

PTMN  (MBOUND.NBPl)  =PFflC-!SUMlNT 

126 

40 

CONTINUE 

127 

DO  45  1=MBP1.M01 

128 

DO  45  J=NBP1,N01 

129 

K»I-1 

130 

L-J-1 

131 

HRITE(6,801)K.L.PTMN(I.J) 

132 

801 

FORMfiTC  P(T,'.  15, '.MS,')  =  '.D17.il) 

133 

45 

CONTINUE 

134 

c 

135 

c 

bflTfi  ADJUSTMENT  FOR  PLOT 

136 

c 

137 

NR0M>41 

138 

NC0L=41 

139 

NKXT=100 

140 

L1NES=0 

141 

flLPHfl=lS. 

142 

BETfl=30. 

143 

FL(1)«0.0 

144 

PL  (2)  =0.0 

145 

SIZE(n=6.0 

146 

SIZE  (2)  =8.0 

147 

c 

148 

c 

SCRLING 

149 

c 

ISO 

PMflX  =  PTMN(l,  1) 
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151  DO  50  1"1.41 

152  Xd)"! 

153  DO  50  J=l,'41 
151*  T(J)=J 

155  IF  (PTMN  (1.  J)  .GT.PMftX)  PMflX  =  PTMN  (I ,  J) 

156  50   CONTINUE 

157  C0NST=8.0/PMflX 

158  DO  51  I>l.m 

159  DO  51  J=1.>11 

160  P(I.J)  =  SNGL  (CONST«PTMN  (I.  J)  ) 

161  51   CONTINUE 

162  C 

163  C 

164  C  IF   OTHER  FORCE  LEVELS  CHANGE  PLOT  flRGUHENTS 

165  C 

166  CALL  PLT301  (X.NROH.T.NCOL.P.flLPHfl.BETfi.FL.TTL.SIZE.WK. 

167  X  ION. KX.KY.NKXT. LINES) 

168  C 

169  C   EXPECTED  VALUES  AND   VARIANCES 

170  C 

171  EM'O.O 

172  EMM=0.0 

173  VARM=0.0 
17q  SPROB-0.0 

175  DO  60  1I=MBP1.M01 

176  1=11-1 

177  SPT=0.0 

178  DO  601  JJ=NBP1,N01 

179  SPT=SPT+   PTMN(ll.JJ) 

180  601   CONTINUE 

181  EM=EM+I«SPT 

182  EMM=EMM+  1«ImSPT 

183  SPR0B=SPRQB*SPT 
1811  60  CONTINUE 

185  VARM»EMM-EM«EM 

186  EN^O.O 

187  ENN=0.0 

188  VARN=0.0 

189  DO  70  II=NBP1.N01 

190  I>II-1 

191  SPT-0.0 

192  DO  701  JJ=MBP1.H01 

193  SPT  =  SPT*   PTMN  (JJ. II) 

194  701   CONTINUE 

195  EN=EN*I«SPT 

196  ENN=ENN+  I«I«SPT 

197  70  CONTINUE 

198  VARN^ENN-ENwEN 

199  WRITE  (6,805)  EH. EN. VRRM. VARN 

200  805  FORMAT ('  ',  'EXPECTED  VALUES  M  .  N".   5QX.'  VARIANCE  M.N  ',/ 
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201  n*    '.On.n.lOX.  017.  11. T51. 017.  11.  lOX.  017.11) 

202  WRITE  (6,1100)  SPROB 

203  1100  FORMflTC  SUM  OF  PROBRBILI T  lES  '.017.11) 

204  STOP 

205  END 
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